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ABSTRACT 

We assess models for the assembly of supermassive black holes (SMBHs) at the center of galaxies that 
trace their hierarchical build-up far up in the dark halo 'merger tree'. Motivated by the recent discovery 
of luminous quasars around redshift z 6 - suggesting a very early assembly epoch - and by numer- 
ical simulations of the fragmentation of primordial molecular clouds in cold dark matter cosmogonies, 
we assume that the first 'seed' black holes (BHs) had intermediate masses and formed in (mini)halos 
collapsing at z ~ 20 from high-cr density fluctuations. As these pregalactic holes become incorporated 
through a series of mergers into larger and larger halos, they sink to the center owing to dynamical 
friction, accrete a fraction of the gas in the merger remnant to become supermassive, form a binary 
system, and eventually coalesce. The merger history of dark matter halos and associated BHs is followed 
by cosmological Monte Carlo realizations of the merger hierarchy from early times until the present in 
a ACDM cosmology. A simple model, where quasar activity is driven by major mergers and SMBHs 
accrete at the Eddington rate a mass that scales with the fifth power of the circular velocity of the host 
halo, is shown to reproduce the observed luminosity function of optically-selected quasars in the redshift 
range 1 < z < 5. A scheme for describing the hardening of a BH binary in a stellar background with 
core formation due to mass ejection is applied, where the stellar cusp oc r~ 2 is promptly regenerated 
after every major merger event, replenishing the mass displaced by the binary. Triple BH interactions 
will inevitably take place at early times if the formation route for the assembly of SMBHs goes back 
to the very first generation of stars, and we follow them in our merger tree. The assumptions underly- 
ing our scenario lead to the prediction of a population of massive BHs wandering in galaxy halos and 
the intergalactic medium at the present epoch, and contributing < 10% to the total BH mass density, 
Psmbh = 4 x 10 5 M Q Mpc~ 3 (h = 0.7). The fraction of binary SMBHs in galaxy nuclei is of order 10% 
today, and it increases with redshift so that almost all massive nuclear BHs at early epochs are in binary 
systems. The fraction of binary quasars (both members brighter than 0.1 L*) instead is less than 0.3% 
at all epochs. The nuclear SMBH occupation fraction is unity (0.6) at the present epoch if the first seed 
BHs were as numerous as the 3.5-cr (4-tr) density peaks at z = 20. 

Subject headings: cosmology: theory - black holes - galaxies: evolution - quasars: general 



1. INTRODUCTION 

Dynamical evidence indicates that supermassive black 
holes (SMBHs) reside at the center of most nearby galax- 
ies (Richstone et al. 1998). The available data show an 
empirical correlation between bulge luminosity and black 
hole mass (Magorrian et al. 1998), which becomes remark- 
ably tight when the stellar velocity dispersion of the host 
bulge, <7 C , is plotted instead of luminosity (Ferrarese & 
Merritt 2000; Gebhardt et al. 2000). The m B H-o- c rela- 
tion implies a rough proportionality between SMBH mass 
and the mass of the baryonic component of the bulge. It 
is not yet understood if this relation was set in primordial 
structures, and consequently how it is maintained through- 
out cosmic time with such a small dispersion, or indeed 
which physical processes established such a correlation in 
the first place (e.g., Silk & Rees 1998; Haehnelt & Kauff- 
mann 2000; Adams, Graff, & Richstone 2001; Burkert & 
Silk 2001). Most recently, it has been shown by Ferrarese 
(2002) that in elliptical and spiral galaxies the bulge ve- 
locity dispersion correlates tightly with the value of the 
circular velocity measured well beyond the optical radius, 



suggesting that ttibh is actually determined by the mass 
of the host dark matter halo. 

The strong link between the masses of SMBHs and the 
gravitational potential wells that host them suggests a fun- 
damental mechanism for assembling black holes and form- 
ing spheroids in galaxy halos. In popular cold dark matter 
(CDM) 'bottom-up' cosmogonies, small-mass subgalactic 
systems form first to merge later into larger and larger 
structures. Galaxy halos then experience multiple mergers 
during their lifetime, with those between comparable-mass 
systems ( "major mergers" ) expected to result in the forma- 
tion of elliptical galaxies (see, e.g., Barnes 1988; Hernquist 
1992). Simple models in which SMBHs are also assumed 
to grow during major mergers and to be present in every 
galaxy at any redshift - while only a fraction of them is 
'active' at any given time - have been shown to explain 
many aspects of the observed evolution of quasars (e.g. 
Cattaneo, Haehnelt, & Rees 1999; Cavaliere & Vittorini 
2000; Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002). 
In hierarchical structure formation scenarios, the ubiquity 
of SMBHs in nearby luminous galaxies can arise even if 
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Assembly and Merger History of SMBHs 



only a small fraction of halos harbor SMBHs at high red- 
shift (Menou, Haiman, & Narayanan 2001). Yet several 
important questions remain unanswered, most notably: 

1. Did the first massive black holes (BHs) form in sub- 
galactic units far up in the merger hierarchy, well before 
the bulk of the stars observed today? The seeds of the 
recently discovery z « 6 quasars in the Sloan Digital Sky 
Survey (SDSS, Fan et al. 2001b) had to appear at very 
high redshift, z > 10, if the SDSS quasars are accreting no 
faster than the Eddington rate and are not gravitationally 
lensed or beamed (Haiman & Loeb 2001). 

2. How massive were the initial BH seeds? A clue to 
this question might lie in the numerous population of ul- 
traluminous off-nuclear ('non-AGN') X-ray sources that 
have been detected in nearby galaxies (e.g. Colbert & 
Mushotzky 1999; Makishima et al. 2000; Kaaret et al. 
2001). Assuming isotropic emission, the inferred masses 
of these 'ULXs' often suggest intermediate-mass BHs with 
m, > a few hundred M Q . 

3. How efficiently do SMBHs and their 'seeds' spiral 
inwards and coalesce as they get incorporated through a 
series of mergers into larger systems? And what is the 
cumulative dynamical effect of multiple BH mergers on 
galaxy cores? 

4. Is there a population of relic 'Population III' massive 
holes lurking in present-day galaxy halos? 

In this paper we explore a formation route for the as- 
sembly of SMBHs in the nuclei of galaxies that traces 
their seeds back to the very first generation of stars, in 
(mini)halos above the cosmological Jeans mass collapsing 
at z ~ 20 from the high-cr peaks of the primordial density 
field. The first stars must have formed out of metal-free 
gas, with the lack of an efficient cooling mechanism possi- 
bly leading to a very top-heavy initial stellar mass function 
(IMF; Larson 1998), and in particular to the production 
of 'very massive stars' (VMSs) with > 100 M Q (Carr, 
Bond, & Arnett 1984). Recent numerical simulations of 
the fragmentation of primordial clouds in standard CDM 
theories all show the formation of Jeans unstable clumps 
with masses exceeding a few hundred solar masses; because 
of the slow subsonic contraction - a regime set up by the 
main gas coolant, molecular hydrogen - further fragmen- 
tation into sub-components is not seen (Bromm, Coppi, & 
Larson 1999, 2002; Abel, Bryan, & Norman 2000). More- 
over, the different conditions of temperature and density 
of the collapsing cloud result in a mass accretion rate over 
the hydrostatic protostellar core ~ 10 3 times larger than 
what observed in local forming stars, suggesting that Pop 
III stars were indeed very massive (Omukai & Nishi 1998; 
Ripamonti et al. 2002). If VMSs form above 260 M , af- 
ter 2 Myr they would collapse to massive BHs containing 
at least half of the initial stellar mass (Fryer, Woosley, & 
Heger 2001), i.e., with masses intermediate between those 
of the stellar and supermassive variety. It has been sug- 
gested by Madau & Rees (2001, hereafter MR) that a nu- 
merous population of massive BHs may have been the end- 
product of the first episode of pregalactic star formation; 
since they form in high-cr rare density peaks, relic massive 
BHs with to. > 150 M Q would be predicted to cluster 
in the cores of more massive halos formed by subsequent 
mergers. 

In this paper we expand upon the original suggestion of 
MR and assess a model for the assembly of SMBHs in the 



nuclei of luminous galaxies out of accreting Pop III seed 
holes. The merger history of dark matter halos and asso- 
ciated black holes is followed through Monte Carlo realiza- 
tions of the merger hierarchy (merger trees). Merger trees 
are a powerful tool for tracking the evolution of SMBH 
binaries along cosmic time, and analyze how their fate is 
influenced by the environment (e.g stellar density cusps). 
We study the conditions under which pregalactic massive 
holes may sink to the halo center owing to dynamical fric- 
tion, accrete a fraction of the gas in the merger remnant to 
become supermassive, form a binary system, and eventu- 
ally coalesce. Major mergers are frequent at early times, so 
a significant number of binary SMBH systems is expected 
to form then. To anticipate the results of our analysis, we 
find that a simple model where quasar activity is driven 
by major mergers and SMBHs accrete at the Eddington 
rate a mass that scales with the fifth power of the cir- 
cular velocity of the host halo, is able to reproduce the 
observed luminosity function of optically-selected quasars 
in the redshift range 1 < z < 5. Minor mergers are largely 
responsible for a population of isolated BHs wandering in 
galaxy halos, while intergalactic BHs will be produced by 
the gravitational slingshot - the ejection of one BH when 
three holes interact. 

2. HALO MERGER TREE 

There are now a number of algorithms for construct- 
ing merger trees, the difficulties and drawbacks of various 
techniques having been reviewed by Somerville & Kolatt 
(1999). We have developed a Monte Carlo algorithm simi- 
lar in spirit to the one described by Cole et al. (2000), and 
based on the extended Press-Schechter formalism (EPS). 
This gives the fraction of mass in a halo of mass M at red- 
shift zo, which at an earlier time was in smaller progenitors 
of mass in the range M to M + dM, 
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where D = S c (z) — 5 c (zq) and S = (j 2 m {z) — <j\ Io {zq) (Bower 
1991; Lacey & Cole 1993). Here <j\j{z) and a 2 M (zo) are 
the linear theory rms density fluctuations smoothed with 
a 'top-hat' filter of mass M and M at redshifts z and z , 
respectively. The 8 c {z) and 5 c (zo) are the critical thresh- 
olds on the linear overdensity for spherical collapse at the 
two redshifts. Integrating this function over the range 
< M < M gives unity: all the mass of M was in 
smaller subclumps at an earlier epoch z > z a . Taking 
the limit z — > zq and multiplying by the factor Mq/M to 
convert from mass weighting to number weighting, equa- 
tion (1) gives the number of progenitors the more massive 
halo fragments into when one takes a small step 5z back 
in time, 
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Our algorithm uses this expression to build a binary 
merger tree that starts from the present day and runs 
backward in time 'disintegrating' a parent halo into its pro- 
genitors. Because for CDM-like power spectra the number 
of halos diverges as the mass goes to zero, it is necessary 
to introduce a cut-off mass or effective mass resolution 
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M rcs , which marks the transition from progenitor - all ha- 
los with M > M ros - to accreted mass - the cumulative 
contribution of all halos with M < M ICS (Somerville & 
Kolatt 1999). Having specified the mass resolution, one 
can compute the mean number of fragments in the range 
M rcs < M < M /2, 



r-Mo/2 m 



and the fraction of accreted mass, 



Jo 



dN M 

dM Mr, 



dM. 
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(4) 



Both of these quantities are proportional to the timestep 
Sz, which is chosen to ensure that multiple fragmentation 
is unlikely, i.e., N p <C 1. Hence it is the frequency of 
mergers that directly controls the timestep: this enables 
the algorithm to follow the merger process with high time 
resolution. Following Cole et al. (2000), at every timestep 
a random number < R < 1 is generated and compared 
to N p . If R > N p the parent halo does not fragment in 
this timestep, but its mass is reduced to account for the 
accreted matter, i.e., a new halo is produced with mass 
Mo(l — -F a )- Fragmentation occurs instead if R < N p : 
then a random value of M in the range M ros < M < Mo /2 
is generated from the distribution in equation (2), to pro- 
duce two new halos with masses M and M (l — F a ) — M. 
The merger hierarchy is built up by repeating the same 
procedure on each subclump at successive timesteps. 

To fully define the tree we need to specify the power- 
spectrum of density fluctuations, which gives the func- 
tion (7m, and the cosmological parameters Slo and Sl\, 
which the critical overdensity for collapse 5 C depends upon. 
Unless otherwise stated, all the results shown below re- 
fer to the currently favoured (by a variety of observa- 
tions) ACDM world model with St = 0.3, Sl A = 0.7, 
h = 0.7, Qbh 2 = 0.02, and n = 1. In this cosmology 
the redshift dependence of the matter density parameter 

and the linear 
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theory growth factor is accurately approximated by 
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(Carroll, Press, & Turner 1992), so that ctm{z) = 
o-m(0)D(z)/D(0). The normalization of the mass fluc- 
tuation spectrum, derived from the abundance of X-ray 
emitting clusters observed in the local universe, is o~s(0) = 
<7 (r = 8ft _1 Mpc) = 0.93 (Eke, Cole, & Frenk 1996). 
We have used the fit to the CDM power-spectrum given 
by Bardeen et al. (1986), modified to account for the ef- 
fects of baryon density following Sugiyama (1995). For the 
spherical collapse density threshold S c we use the fit (ac- 
curate within 0.1% for 1 — SI > 0.01) from Nakamura & 
Suto (1997), S c (z) = 1.686[1 + 0.012299 log(l - SI)]. 

All that remains to be fixed is the mass resolution and 
timestep. We have taken M res to represent at z — a 
subclump with virial velocity equal to 10% that of the Mq 
parent halo, i.e., M res = 10~ 3 M : the mass resolution is 
taken to decrease with redshift as (1 + z)~ 3 - 5 , so it is al- 
ways less than 5% of the mass of the main halo in the 



merger hierarchy. This ensures a sufficiently wide range of 
masses to allow both for minor and major mergers in the 
tree at all redshifts. Due to the binary nature of the algo- 
rithm (which allows only a single fragmentation within a 
timestep), M res cannot be chosen to be arbitrarily small, 
as in this case multiple fragmentation becomes unprobable 
(i.e., N p <C 1) only if the time resolution (and hence the 
computational time) is extremely high. For instance, if we 
keep the ratio between M res and the mass of the main halo 
in the tree fixed at the present-day value, the time for one 
Monte Carlo realization is longer by an order of magni- 
tude. Since our aim is to track the merger hierarchy to 
very high redshifts, to keep the computational time down 
to acceptable values, we have approximated the fraction 
of accreted mass by the fitting formula 



F a = a[log(M/M r( 



(6) 



instead of integrating numerically the mass function of 
progenitors with M < M res at every branch of the tree. 
Here, depending on mass and redshift, a and b are in the 
range (0.3 — 40) x 10~ 3 and 0.3 — 0.8, respectively, and have 
been adjusted empirically to provide a fit to equation (4) 
within 5% in the range of masses and reshifts considered. 
Beside conserving mass, a merger tree algorithm must re- 
produce at all redshifts the conditional mass function pre- 
dicted by the EPS theory The comparison is shown in 
Figures 1 and 2: for this set of realizations we have used 
820 timesteps logarithmically spaced in expansion factor 
between z = and z = 20. With the above prescriptions, 
the tree typically agrees with the EPS predictions within 
a factor of 2 up to z = 20 for masses greater than 3xM res . 
Each realization we generate tracks backwards the merger 
history of 220 parent halos at the present epoch picked 
over the mass range 10 11 < M < 10 15 M Q , where the 
lower limit has been chosen to match the minimum mass 
effectively probed by dynamical studies of SMBH hosts in 
the local universe. 




Fig. 1. — Mean number of progenitors with mass M for a z$ = 0, 
Mo = 10 11 Mq parent halo, at redshifts z = 1,5, 15, 20. Solid lines: 
predictions of the EPS theory. Histograms: results for the merger 
tree (mean of 50 realizations), M > 2 X M TCB . Error bars represent 
the Poissonian error in the counts. 
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These are broken up into as many as 70,000 progeni- 
tors by z = 20. One issue of concerns involves systematic 
deviations of the unconditional (and conditional) Press- 
Schechter (PS) mass function compared with N-body sim- 
ulations. At low rcdshifts, the PS mass function overpre- 
dicts the number of small halos by a factor of 1.5 to 2 (e.g. 
Gross et al. 1998; Sheth & Tormcn 1999). The model and 
simulation results agree well on all scales at z <~ 1, while at 
higher rcdshifts the abundance of large-mass halos is un- 
derestimated by the PS model (Somerville et al. 2000). In 
general, halo merging histories constructed using the EPS 
formalism are reasonably consistent with those extracted 
from N-body simulations. At very high redshift (z > 10), 
a recent cosmological simulation by Jang-Condcll & Hern- 
quist (2001) finds good agreement with the PS mass func- 
tion for the lowest halo masses of interest here. 
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Same as Fig. 1 for a z = 0, M = 4 X 10 13 M parent 



3. GROWTH OF SUPERMASSIVE BLACK HOLES 
3.1. Pregalactic 'seeds' 

Following MR, we will assume that one seed BH forms 
in each of the rare density peaks above v-a at z = 20. 
All seed BHs are assigned a mass m, = 150 M Q . In gen- 
eral, we find that our results are not very sensitive to the 
precise mass choice for these seed holes (see §4). As our 
fiducial model we take v = 3.5, corresponding in the as- 
sumed ACDM cosmology to minihalos of mass M secd = 
1.1 x 10 7 hr 1 M©. This is larger than the minimum mass 
threshold for baryonic condensation, M m i n «5x 10 5 M 
(ft = 0.65), found in the numerical simulations of Fuller 
& Couchman (2000). Above M m i n the H 2 cooling time 
is shorter than the Hubble time at virialization, the gas 
in the central halo regions becomes self-gravitating, and 
stars can form. We notice that the EPS formalism pre- 
dicts of order 1 progenitor above M see d for our lower mass 
parent halo, M = 10 11 M©. Also, our resolution mass, 
M res = 10~ 3 M (1 + z)- 3 ' 5 is always lower at z = 20 than 
M S eed for all parent halos with masses Mo < 10 15 M©. 

A pregalactic halo at z — 20 is characterized by a virial 
radius (defined as the radius of the sphere encompassing 

4 For the assumed cosmology this can be approximated by A v ; r = 



a mean mass density A v i r /5 cr it, where p cr ;t is the criti- 
cal density for closure at the redshift z and A v j r is the 

density contrast at virialization 4 ) r v - lr = 390 pc ft, -1 , 

and a circular velocity V c — 10.5 kms -1 M^ 3 at r v j r . 
The gas collapsing along with the dark matter pertur- 
bation will be shock heated to the virial temperature 



2/3 

3200 K M 7 , where is the halo mass in units 
of 10 7 ft _1 M©. The total baryonic mass within the virial 
radius is equal to (M see d^fc/^o)- In a Gaussian theory, 
halos more massive than the v-a peaks contain a fraction 
erfc(z//\/2) (= 0.00047 for v = 3.5) of the mass of the uni- 
verse. Therefore the mass density parameter of our '3.5-er' 
pregalactic holes is 



n. = 



0.00047 fi m. 
1.1 x 10 7 ft-! M© 



> 2 x 10" 9 ft. 



(7) 



This is much smaller than the density parameter of the 
supermassive variety found in the nuclei of most nearby 
galaxies, 

_ 6 x 10 5 ft M© Mpc- 3 



'SMBH = 



" 1V1 Q 1V1 ^ «2xl0- 6 ft- 1 (8) 

2.8 x IO^MqMpc- 3 ' 1 ' 



where the value at the denominator is the critical density 
today and we have taken for the local density of SMBHs 
the value recently inferred by Merritt & Ferrarese (2001). 
Clearly, if SMBHs form out of very rare Pop III BHs, the 
present-day mass density of SMBHs must have been accu- 
mulated during cosmic history via gas accretion, with BH- 
BH mergers playing a secondary role. This is increasingly 
less true, of course, if the seed holes are more numerous 
and populate the 2- or 3-cr peaks instead, or halos with 
smaller masses at z > 20 (MR). 

The choice of where to initially locate our seed BHs, 
while motivated by recent numerical simulations of the 
formation of the first stars, is clearly somewhat arbi- 
trary. Computational costs do not allow us to follow 
the merger hierarchy much beyond z = 20, or down to 
minihalo masses smaller than M scc( j. As argued by MR, 
if an extreme IMF is linked to primordial H 2 chemistry 
and cooling, it seems unlikely that the formation of mas- 
sive BHs from zero-metallicity VMSs might have been a 
very efficient process, since metal-free VMSs are copious 
sources of Lyman- Werner photons (Bromm, Kudritzki, & 
Loeb 2001). This radiation may photodissociate H 2 else- 
where within the host halo, escape into the intergalactic 
medium, and form a soft UV background that suppresses 
molecular cooling throughout the universe (Haiman, Abel, 
& Rees 2000) and may inhibit the formation of a much 
more numerous population of VMSs. Moreover, VMSs in 
the mass range 140 < < 260 M© are predicted to 
make pair-instability supernovae with explosion (kinetic) 
energies of up to 10 53 erg. This is typically much larger 
than the baryon binding energy of a subgalactic fragment, 
E b w 10 52 - 5 erg h^flbM^ 3 [(1 + z)/20]. Minihalos will 
then be completely disrupted ('blown-away') by such en- 
ergetic events, metal-enriched material will be returned to 
the IGM (e.g. Madau, Ferrara, & Rees 2001) and later 
collect in the cores of more massive halos formed by sub- 
sequent mergers, where a 'second' generation of stars may 

178Q - 45 (Eke, Navarro, & Frcnk 1998). 
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now be able to form with an IMF that is less biased to- 
wards very high stellar masses (Schneider et al. 2002). 

3.2. Major mergers 

We model each dark halo as a singular isothermal sphere 
(SIS) with circular velocity V c , one-dimensional velocity 
dispersion <tdm = V c /V2, and density p(r) = V c 2 /4irGr 2 , 
truncated at the virial radius. When two halos of mass 
M and M s merge, the 'satellite' (less massive) progen- 
itor (mass M s ) is assumed to sink to the center of the 
more massive pre-existing system on the Chandrasekhar 
dynamical friction (against the dark matter background) 
timescale 

r 2 . V 1 4- P 1 

tdf = 1-17 „ " rc , \ £ ° = 1-65—^- — tt — -e (9) 
GM S In A P Hy/K~\nA y ' 

(Lacey & Cole 1993; Binney & Tremaine 1987), where V c 
is the circular velocity of the satellite in the new halo of 
mass M+M s and virial radius r vir , r circ is the radius of the 
circular orbit having the same energy as the actual orbit, 
the 'circularity' e is the ratio between the orbital angu- 
lar momentum and that of the circular orbit having the 
same energy, H is the Hubble parameter, P = M s /M is 
the (total) mass ratio of the progenitors, and the Coulomb 
logarithm is taken to be In A m ln(l + P). The dependence 
of this timescale on the orbital parameters is contained in 
the term 

e = e"(r circ /r vir ) 2 . (10) 

The most likely orbits occurring in cosmological CDM sim- 
ulations of structure formation have circularity e = 0.5 
and r c i rc /r v ; r = 0.6 (e.g. Tormcn 1997; Ghigna et al. 
1998). With these initial orbital parameters, recent nu- 
merical investigations by van den Bosch et al. (1999) and 
Colpi, Mayer, & Governato (1999) suggest a value a = 0.4- 
0.5 for the exponent in equation (10). Here we assume 
= 0.3, but we note that the merger timescale computed 
in this way does not include the increase in the orbital de- 
cay timescale due to tidal stripping of the satellite (Colpi 
et al. 1999). Satellites will merge with the central galaxy 
on timescales shorter than the then Hubble time only in 
the case of major mergers, P > 0.3. In minor mergers 
tidal stripping may leave the satellite BH almost 'naked' 
of its dark halo, too far from the center of the remnant for 
the formation of a black hole binary. 

Figure 3 shows the number of major mergers per unit 
rcdshift bin experienced by halos of different masses. For 
galaxy-sized halo this quantity happens to peak in the red- 
shift range 2-4, the epoch when the observed space density 
of optically-selected quasar also reaches a maximum. Hy- 
drodynamic simulations of major mergers have shown that 
a significant fraction of the gas in interacting galaxies falls 
to the center of the merged system (Mihos & Hernquist 
1994, 1996): the cold gas may be eventually driven into 
the very inner regions, fueling an accretion episode and the 
growth of the nuclear BH. In the following we shall make 
the simplifying assumption that SMBHs accrete material 
only during major mergers. 
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Fig. 3. — Mean number of major mergers experienced per unit 
rcdshift by halos with masses > 10 10 Mq. Solid line: progenitors 
of a Mo = 2 X 10 13 Mq halo at z = 0. Dashed line: same for 
M = 3 X 10 12 Mq. Dotted line: same for M = 10 11 Mq. 

3.3. Accretion history 

The physical processes that determine the amount of ac- 
creted gas and the characteristic accretion timescales onto 
SMBHs are poorly understood, and different prescriptions 
have been proposed in the literature to explain the ob- 
served evolution of QSOs within hierarchical clustering 
cosmologies. We shall not attempt here to model these 
processes in details, but we note that the fraction of cold 
gas ending up in the hole must depend on the properties 
of the host halo in such a way to ultimately lead to the ob- 
served correlation between stellar velocity dispersion and 
SMBH mass. Using the most up-to-date set of black hole 
mass measurements, Ferrarese (2002) finds 

m BH = (4.4 ± 0.9) x 10 7 M <f 5 8 ±0 - 52 , (11) 

where <7 Cl i5n is the bulge velocity dispersion (defined within 
an aperture of size < 0.5kpc) in units of 150 kms -1 . Geb- 
hardt et al. (2000) and Tremaine et al. (2002) report a sim- 
ilar relation with a somewhat shallower slope. From a sam- 
ple of spirals and elliptical galaxies with cr c > 70 kms -1 , 
Ferrarese (2002) also shows that the stellar velocity dis- 
persion is strictly correlated with the asymptotic value of 
the circular velocity V c measured well beyond the optical 
radius, 

log V c = (0.88 ±0.17) log cr c + (0.47 ±0.35). (12) 

To avoid introducing additional parameters to our model, 
as well as uncertainties linked to gas cooling, star forma- 
tion, and supernova feedback, we combine the two previous 
relations and adopt the following simple prescription for 
the mass accreted by a SMBH during each major merger: 

Am acc = 1.3 x 10 4 M Q K V c % , (13) 

where V Ci iso is the circular velocity of the merged system 
in units of 150 kms -1 (cf. Kauffmann & Hachnelt 2000). 
We assume this mass is accreted by the BH in the more 
massive progenitor halo at the Eddington rate, after about 
a dynamical timescale (estimated at one tenth the virial 
radius, td yn — OAr^/Vc). The normalization factor K 
is of order unity and is fixed in order to reproduce the 
TObh — cr c relation observed locally. 
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The above relations are only valid up to a value of the ve- 
locity dispersion corresponding to galaxy group and clus- 
ter scales. For instance, the SMBH in M87 has a mass of 
TObh = 3 x 10 9 M Q (Harms et al. 1994), perfectly corre- 
lating with the observed circular velocity V c = 506 km s -1 
of the galaxy, but not with the circular velocity of the 
Virgo cluster (V c w lOOOkms -1 ). Moreover, the larger 
halo (V c = 1280 km s -1 ) considered in our merger-tree 
set would contain at z — a SMBH more massive than 
10 10 M Q : this would lead to an overestimate of the quasar 
LF at late epochs. Following Kauffmann & Haehnclt 
(2000), we have then inhibited gas accretion onto SMBHs 
in all halos with V c > 600 km s -1 . We find that this as- 
sumption affects only the accretion history of the SMBHs 
hosted in the two more massive halos of our realizations, 
and only at late epochs. On small mass scales, on the 
other hand, accretion onto BHs hosted by minihalos with 
virial temperature T vlr < 10 4 K may be inhibited by UV 
radiation. Photoionization by an internal UV source or a 
nearby external one will heat the gas inside these shallow 
potential wells to 10 4 K, leading to the photoevaporation 
of baryons out of their hosts. After the reionization epoch 
gas cooling, star formation, and accretion onto BHs may 
only possible within more massive halos with virial tem- 
peratures T v i r > a few x 10 4 K, where pressure support is 
reduced and gas can condense due to atomic line processes. 
In our model we therefore suppress gas accretion onto all 
BHs in minihalos with T v ; r < 10 4 . We have checked that 
this assumption has little or no effect at z <J 10. 

Modeling gas accretion onto BHs with the recipes just 
described we find that, along cosmic history, most of the 
final mass of SMBHs come from gas accretion, rather than 
from BH merging. The final mass m of the SMBH formed 
after coalescence assumes the entropy-area relation for 
BHs (maximally efficient radiative merging, e.g. Ciotti 
& van Albada 2001): the total entropy S of the system 
remains unchanged, S = m 2 /4 = Si + S2 = m 2 /4 + m|/4 
(taking G = c = k = h = 1; Hawking & Ellis 1973). In 
contrast with previous work (e.g. Kauffmann & Hachnelt 
2000; Menou et al. 2001), we do not make here the simpli- 
fying assumption that the two pre-existing holes coalesce 
instantaneously. The evolution of SMBH pairs will be dis- 
cussed in the next section. 



3.4. Dynamical evolution of BH binaries 

In our model the merging - driven by dynamical friction 
against the dark matter background - of two halo+BH sys- 
tems with mass ratio P > 0.3 will drag in a satellite BH 
towards the center of the more massive progenitor; this will 
inevitably lead to the formation of a bound SMBH binary 
in the violently relaxed core of the newly merged stellar 
system. Figure 4 shows the typical mass ratio of binary 
BHs together with the mean mass of the larger member of 
the pair as a function of redshift. At late epochs most of 
the BH pairs have unequal masses and, since the growth 
history of SMBHs does not track that of dark matter halos, 
a major merger between halos does not necessarily result 
in a BH binary with a large P, as evident from the upper 
panel of the figure. 




Fig. 4. — The mass ratio of BH binaries (upper panel) and the 
mass of the larger BH (lower panel) are shown as a function of red- 
shift. The points show the mean values and 1-<t error bars from all 
Monte Carlo realizations. At high redshift binary members are seed 
Pop III BHs with nearly equal masses; as time goes on the holes 
grow mainly due to gas accretion, and low mass ratios become more 
probable. 

The subsequent evolution of BH binaries was first out- 
lined by Begelman, Blandford, & Rees (1980). Consider a 
binary with BH masses mi > m 2 and semimajor axis a(t) 
in an isotropic background of stars of mass m* <C m 2 and 
density p±(r). We use a simple model for the initial cen- 
tral stellar distribution, an SIS with a velocity dispersion 
comparable to the halo itdm, 



"DM 

2irGr 2 ■ 



(14) 



This appears to be a good assumption for early-type lens 
galaxies (e.g. Koopmans & Treu 2002). When the age of 
the system is larger than the stellar relaxation time, the 
equilibrium distribution of stars around a BH is expected 
to be cuspy, oc r~ 7 / 4 , within the gravitational sphere of 
influence of the BH, even if the original profile had a core 
(Bahcall & Wolf 1976). 

The binary forms at a separation = G(mi + 
ra2)/(2(TQ M ) at which the enclosed stellar mass equals 
mi + m2, and initially hardens by dynamical friction 
from distant stars acting on each BH individually. But 
as the binary separation shrinks (the binary 'hardens') 
the effectiveness of dynamical friction slowly declines be- 
cause distant encounters perturb only the binary center's 
of mass but not its semimajor axis. The BH pair then 
hardens via three-body interactions, i.e., by capturing the 
stars that pass within a distance ~ a of it and ejecting 
them ('gravitational slingshot') at much higher velocities, 
v e j ~ Vbin = [G(mi + m^/ct] 1 / 2 , where Vt>i n is the relative 
velocity of the two BHs if their orbit is circular: this is the 
hard binary stage. In Quinlan's (1996) simulations of the 
dynamical evolution of massive BH binaries, the system 
does not become hard until a falls below 



ah = 



Gm 2 



1 pc 



m 2 



10 7 - 3 M Q 



-2 

T DM.150 



(15) 
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(Quinlan 1996). 5 We assume that the 'bottleneck' stages 
of the binary shrinking occur for separations a < a^; in 
a major merger, after a dynamical friction timescale, we 
form the BH binary at a separation ah and let it evolve. 

In a fixed background, the hardening timescale \a/a\ de- 
creases with a, 



t h = 



G Pi ,aH 



= 10 34 yr 



pc 



'DM.150' 



(16) 



the problems associated with loss-cone depletion (which 
may ultimately cause the binary to 'stall') and helps the 
binary merge (Quinlan & Hcrnquist 1997; Milosavljevic & 
Merritt 2001). Here we adopt a simple analytical scheme 
following Merritt (2000) that qualitatively reproduces the 
evolution observed in N-body simulations. If M. c ) is the 
stellar mass ejected by the BH pair, the binary evolution 
and its effect on the galaxy core are determined by the 
coupled equations 



and the binary would spend the longest period of time 
with a Ki ah- Here the second equality assumes an SIS (eq. 
14) down to a distance a from the center, and the dimen- 
sionlcss hardening rate is H k, 15 in the limit of a very 
hard, equal-mass binary (Quinlan 1996). If the harden- 
ing continues sufficiently far, gravitational radiation losses 
can take over, and the two BHs rapidly coalesce on the 
timescale (for a circular orbit) 



t _ 5c 5 a 4 (t) 

gr 256G 3 TOiTO2(mi + rri2) 



(17) 



(Peters 1964). If the binary can shrink to a separation 



0.014pc 



(mi + m 2 )mim 2 



10 21 - 3 UJ 



1/4 



(18) 



the binary will coalesce within 10 Gyr due to the emis- 
sion of gravitational waves. Here we have normalized 
to the case mi = m 2 — 10 7 M . In our model we ne- 
glect the recoil due to the non-zero net linear momentum 
carried away by gravitational waves in the coalescence of 
two unequal mass BHs ('gravitational rocket'). Radiation 
recoil is a strong-field gravitational effect that depends 
on the lack of symmetry in the system, and may eject 
massive BHs from galaxy cores. To date, the outcome 
of a gravitational rocket remains uncertain. Newtonian- 
approximation and perturbative calculations of two orbit- 
ing Schwarzschild holes suggest maximum recoil velocities 
of order 70 - 100 kms" 1 (Fitchett 1983, Fitchett & De- 
tweiler 1984). 

In practice, however, one cannot assume a fixed stellar 
background in estimating the rate of BH mergers, as the 
hardening of the binary modifies the stellar density p* in 
equation (16): the shrinking of the pair removes mass inte- 
rior to the binary orbit, depleting the galaxy core of stars 
and slowing down further hardening. The effect of loss- 
cone depletion (the depletion of low-angular momentum 
stars that get close enough to extract energy from a hard 
binary) is one of the major uncertainties in computing 
the merger time, and makes it difficult to construct viable 
merger scenarios for BH binaries. A recent analysis by Yu 
(2002) has shown that in significantly flattened or triaxial 
galaxies the supply of low-angular momentum stars may 
be sufficient to reach a gr . A massive gaseous disk sur- 
rounding the binary may further speed up the merger rate 
(Gould & Rix 2000). N-body simulations of BH binary de- 
cay suggest that the wandering of the binary center of mass 
from the galaxy center (and to a lesser extent the diffu- 
sion of stars into the loss cone) may also work to mitigate 



dfl\ =H G BL 
dt \aj (Tdm 



and 



dM c 



dln(l/a) 



J {mi + m 2 ), 



(19) 



(20) 



where J is the dimcnsionless mass ejection rate, J rj 1 
nearly independent of a for a <C ah (Quinlan 1996). In- 
tegrating the second equation one finds M c j w J(mi + 
m 2 ) \n(ah/a): the binary ejects of order its own mass in 
shrinking from ah to ah/3. We assume that the stellar 
mass removal creates a core of radius r c and constant den- 
sity p c = p*{r c ), so that the ejected mass can be written 
as 



G 



4ir 

-JPc 



4 <T DM r c 

3 G 



(21) 



(Merritt 2000). From equations (20) and (21) one derives 



r c (t) 



4ff DM 



GJ(mi + m 2 ) \n(a h /a), 



and the core density decreases as 



Pc(t) 



8^DM 



9nG 3 M 2 ci (t)' 



(22) 



(23) 



The above relations, assuming a constant ctdm during the 
hardening of the binary, are strictly valid only if the stel- 
lar relaxation timescale is long compared to the hardening 
time. The binary separation quickly falls below r c and 
subsequent evolution is slowed down due to the declining 
stellar density, with a hardening time, 



t h = 



2?rr 2 



(24) 



which now becomes increasingly long as the binary shrinks. 
In this model the mass ejected increases logarithmically 
with time, and the binary can 'heat' background stars at 
radii r, a <C r < r c . In N-body simulations this may hap- 
pen due to the Brownian motion of the binary induced 
by continuous interactions with other stars. Also, stars 
on eccentric orbits are most likely to interact with the bi- 
nary and be removed, then loose their kinetic energy to 
the background as they spiral back in and are kicked out 
again. 

Figure 5 shows the evolution of a binary of seed, 
intermediate-mass BHs in two dark matter halos with 
different velocity dispersion. The binary separation, a, 



5 The standard definition of a "hard" binary, one where its binding energy E\> = Gmim^la exceeds the typical kinetic energy of the sur- 
rounding stars 3m*<TQ M /2 (Binney & Tremaine 1987), is inapplicable to massive BH binaries as they are always hard if bound. Quinlan (1996) 
defines hardness instead in terms of the binary orbital velocity; in his definition a hard binary hardens at a constant rate. 
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shrinks as the pair interacts with the surrounding stellar 
field, and, at the same time, the ejection of stars decreases 
the central density creating a stellar core. 
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Fig. 5. — The evolution (from bottom to top) of core radius, 
core stellar density, and hardening timescale during the shrinking 
of a 'Pop IIP BH binary against its separation. As the pair of 
mm = 150 M0 BHs shrinks, the initial p* oc r~ 2 stellar cusp is 
gradually converted into a constant density core by the gravita- 
tional slingshot, and the hardening timescale lengthens. Thin lines: 
udm = 7 kms -1 . The total stellar mass ejected prior to coalescence 
is 7.5(mi +012). Thick lines: <tdm = 20 kms -1 . The total stellar 
mass ejected prior to coalescence is 6.5(mi + 7112). 

Finally, as it is conceivable that major mergers be- 
tween galaxies may trigger bursts of star formation (e.g. 
Somcrvillc, Primack, & Faber 2001), we further assume 
that a stellar cusp oc r~ 2 is promptly regenerated after ev- 
ery major merger event, replenishing the mass displaced 
by the BH binary. For a fixed binary mass the coa- 
lescence timescale is shorter in the case of more mas- 
sive galaxies. The evolution of two SMBH binaries in a 
(7dm = 200 kms -1 halo is depicted in Figure 6. 

An equal mass binary with m\ — m 2 = 10 7 M Q needs a 
longer time to merge than a binary with mi = 10 7 M Q 3> 
m 2 , as it must eject a larger number of stars. A compar- 
ison with a straightforward extrapolation of Milosavljevic 
& Merritt (2001) N-body results shows that the scheme we 
adopt tends to overestimate the binary evolution timescale 
by about a factor 3 (scaling to real galaxies like M32 and 
M87). Two main factors contribute to this discrepancy: 
first, Milosavljevic & Merritt (2001) let the slope of the 
density profile change smoothly during the hardening of 
the binary, ending with a shallow cusp oc r -1 rather than 
with a flat core; and second, they take into account the 
Brownian motion of the binary, which makes the BHs in- 
teract with a larger number of stars in the central region. 
Nevertheless, the coalescence timescales estimated in our 
model may be too short, since it is the total stellar den- 
sity that is allowed to decrease following equation (23), 
not the density of low angular-momentum stars (i.e. we 
neglect the depopulation of the loss-cone); furthermore 
our model assumes replenishment of the stellar cusp af- 
ter every major merger. Yu (2002) has recently studied 
the merger of binary SMBHs assuming the central stellar 



profiles observed by Faber et al. (1997) in a sample of 
local galaxies, and finds that BH coalescence timescales 
may in some case be longer than the Hubble time. 
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Fig. 6. — The evolution (from bottom to top) of core radius, 
core stellar density, and hardening timescale during the shrinking 
of a SMBH binary in a halo with ctdm = 140 kms -1 . Thin lines: 
mi = m.2 = 10 7 M . The total stellar mass ejected prior to coa- 
lescence is 5(mi + 1712)- Thick lines: mi = 10 7 Mq, W2 = 10 s Mq. 
The total stellar mass ejected prior to coalescence is 2(mi + m^). 



3.5. Triple BH interactions 

The dynamical evolution of SMBH binaries may be dis- 
turbed by a third incoming BH, if another major merger 
takes place before the pre-existing binary has had time 
to coalesce (e.g. Hut & Rees 1992; Xu & Ostriker 1994). 
In a minor merger the intruder BH is stripped of most of 
the surrounding dark and luminous matter; the ensuing 
long dynamical friction timescale does not allow a close 
encounter between the central binary and the intruder. 
Within our scheme these BHs remain wandering in galaxy 
halos through successive mergers. If the incoming hole 
reaches the sphere of influence (determined in our model 
by the hardening distance ah) of the central binary, the 
three BHs are likely to undergo a complicated resonance 
scattering interaction, leading to the final expulsion of one 
of the three bodies (gravitational slingshot). Typically 
an encounter between an intruder of mass m^t smaller 
than both binary members leads to a scattering event, 
where the binary recoils by momentum conservation and 
the incoming lighter BH is ejected from the galaxy nu- 
cleus. The binary also becomes more tightly bound, each 
such encounter typically increasing its binding energy 
by the amount (AE/Eb) ~ 0.4mi nt /(mi + m 2 ) (Hills & 
Fullerton 1980; Colpi, Possenti, & Gualandris 2002). By 
contrast, when the intruder is more massive than one or 
both binary components, the probability of an exchange 
is extremely high: the incoming hole becomes the mem- 
ber of a new binary, and the lightest BH of the original 
pair gets ejected (Hills & Fullerton 1980). For mass ratios 
mint /(mi + m 2 ) <; 2, most of the increase in the binding 
energy of the pair is due to the actual shrinking of the or- 
bit, while above this value the binding energy rises mainly 
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due to the replacement of a low mass member by a more 
massive BH. In the latter case, and for head-on collisions 
and equal mass binaries, the fractional increase in binding 
energy is approximately constant with a value of 3.1 (Hills 
& Fullerton 1980). If the binary is hard the kinetic en- 
ergy and momentum of the intruder are much lower than 
the orbital binding energy and the recoil momentum of 
the ejected body (we have checked a posteriori that this 
is true in nearly all cases of triple interactions) . Conser- 
vation of energy and momentum in the interaction allows 
then to estimate the recoil velocity of the binary and in- 
truder. Let 772 c j be the mass of the lightest of the three 
BHs, and m^in the mass of the final binary, i.e m c j = m^t 
and mbin = m i + m 2 i n a scattering event, m C j = m 2 and 
mbin = mint + mi for exchanges. The kinetic energy of 
the ejected BH and of the binary after the encounter will 
then be 



-Kbin = 



AE 



1 + (mbin/m e j) 



(25) 



K , = ^ . 

CJ l + (m c j/m b in) 



(26) 
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Fig. 7. — The number of triple interactions for different ratios 
between the intruder BH mass m; nt and the mass m\ +n%2 of the bi- 
nary, in different redshift intervals. Left panel: ctdm = 250 kms -1 . 
Right panel: ctdm = 150 kms -1 . This histogram includes results 
from all 20 Monte Carlo realizations of the same halo mass. At 
very high redshift {upper panel, 7 < z < 20) equal mass system 
interactions are more common, while at low redshift (lower panel, 
< z < 3) a high-mass binary typically interacts with an intruder 
of much lower mass. At intermediate-high redshift (middle panel, 
3 < z < 7) a transition regime occurs. 



We adopt a simple scheme, where: 

• if mi n t < m 2 a scattering event occurs, with 
(AE/E b ) = 0.4m int /(mi + 777 2 ). The new semi- 
major axis is ai = a /[l + 0.4 m; n t/(mi + m 2 )]. 



• if m 2 < mint < 2 (toi + 777 2 ) an exchange with 
(AE/Eb) = 0.4mi nt /(mi + m 2 ) takes place. The 
new semi-major axis is a\ = ao(mi n t/m 2 )[l + 

0.4 777 in t/(mi + m 2 )]. 



• if mi n t > 2 (toi + m 2 ) again an exchange hap- 
pens, with (AE/Eb) = 0.9 (we have rescaled the 
value of 3.1 given by Hills & Fullerton 1980 for 
head-on collisions to account for a distribution of 
impact parameters). The new semi- major axis is 
ai = 0.53ao(m in t/m 2 ). 



In all cases we have used the equations (25) and (26) 
to estimate the recoil velocities. All triple interactions are 
followed along the merger tree, as they modify the binary 
separation during each encounter. At high redshift we find 
that the increase in binding energy causes the binary to 
shrink to a separation small enough that coalescence by 
gravitational radiation occurs, since most of these encoun- 
ters happen among approximately equal-mass systems. At 
later epochs events with low m; n t/(mi + m 2 ) mass ratios 
are more common (see Figure 7), and the binding energy 
increases only slightly after such interactions. 



What happens to a BH pair+intruder system after the 
slingshot mechanism? In a SIS+core halo, the gravita- 
tional potential is 



$(r) = -2cr 
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r > r v 



(27) 

where the core radius r c is that created by the hardening 
of the binary at the time of the triple interaction. If the 
kick velocity of the binary and/or single BH exceeds the 
escape speed v csc — y/2\4>\, the hole(s) will leave the galaxy 
altogether. We find that the recoil velocity of the single 
hole is larger than u osc in 99% of encounters. The binary 
is ejected instead in only 8% of the encounters (Figure 8) 
and typically at very high redshifts, when all BHs are in 
the same mass range. 

For equal mass holes both the binary and the single BH 
are ejected from radius r- m to infinity when the orbital ve- 
locity Vbin of the binary satisfies the condition 



Vbin > 7.7v^( 



(28) 



The 1% of single BHs not escaping their host ha- 
los are typically slung to the periphery of the galaxy 
with consequently long dynamical friction timescales; 
most of the binaries recoil instead within the core and 
fall back to the center soon afterwards, with tdf < 
0.01 Gyr. Since most of the ejected holes escape their 
hosts, in our scheme the number of wandering BHs 
due to the gravitational slingshot and retained within 
galaxy halos is thus significantly lower (by about a 
factor of 50) than that left over by minor mergers. 
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Fig. 8. — The fraction of BH binaries (upper panel) and single BHs 
(lower panel) with recoil velocity v larger than the escape speed, as a 
function of redshift. The histogram includes results from all Monte 
Carlo realizations. 



4. IMPLICATIONS 

In this section wc discuss some of the consequences pre- 
dicted by our fiducial scenario (and a few variants) for the 
growth of SMBHs in the nuclei of galaxies. 



4.1. The quasars luminosity function 

In our framework quasar activity is triggered by major 
mergers and SMBHs accrete at the Eddington rate, ttie = 
47rGra p mBH/(c0'Te), where e is the radiation efficiency. 
Accretion starts after about one dynamical timescale and 
lasts until a mass given by equation (13) has been added 
to the hole. Rest mass is converted to radiation with a 
10% efficiency; only a fraction ,/b = 0.08 of the bolomeric 
power is radiated in the blue band. We have compared 
theoretical luminosity functions (LF) at different redshifts 
with the most recent determination of the quasar blue LF 
from the 2dF survey (0.3 < z < 2.3, Boyle et al. 2000) and 
the SDSS (3.3 < z < 5, Fan et al. 2001a). The 2dF LF is 
a double power-law, which we have extrapolated beyond 
redshift 2 assuming pure luminosity evolution; the best fit- 
ting parameters for a ACDM cosmology are given by Boyle 
et al. (2000). The SDSS samples only the very bright end 
of the LF; Fan et al. (2001a) fit a single power-law to the 
data. 

A detailed comparison at very early times required the 
use of ad hoc merger trees to simulate very luminous 
quasars in very massive, rare dark halos at high redshifts. 
For this purpose we have simulated the merging histories 
of 10 15 M Q parent halos starting from z — 3.5, and ap- 
plied the evolutionary scheme for their SMBHs outlined 
in the previous sections. As shown in Figure 9, our sim- 
ple model reproduces reasonably well the observed LF of 
optically-selected quasars in the redshift range 1 < z < 5. 
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Fig. 9. — The B-band luminosity function of quasars at vari- 
ous redshifts. Filled circles: the values predicted from our SMBH 
assembly and merging model history. Error bars indicate the pois- 
sonian error in the counts. Solid lines at z = 1,2: 2dF LF. Solid 
lines at z = 3,5: SDSS LF. The dashed lines show the extrapola- 
tion to faint magnitudes of the best fit LF from Boyle et al. (2000) 
(z = 1, 2), and Fan et al. (2001a) (z = 3, 5). 

The slope at low luminosities matches the one inferred 
by Boyle et al. , and is considerably flatter than the ex- 
trapolation of the SDSS power-law. The tendency to over- 
estimate the number of bright QSOs at z < 1 may be less 
severe due to the presence of a substantial population of 
(optically) obscured luminous AGNs at low redshifts sug- 
gested by recent Chandra results (e.g. Rosati et al. 2002; 
Barger et al. 2001), which may help reducing the discrep- 
ancy. 

To assess the impact on our results of making seed BHs 
more common or rarer, we have run realizations that place 
m, = 150 M Q BHs at z = 20 in 3-cr (lower bias) and 4- 
o~ (higher bias) peaks instead of the fiducial 3.5-cr. We 
find that we can still reproduce the observed quasar LF 
luminosity function by changing the major merger thresh- 
old, from P > 0.3 in the lower bias case, to P > 0.1 in 
the higher bias case. We have also run a case where the 
initial, seed BH mass is m. = 1000 M Q instead of 150 
M©. We find little change at z < 5. The number of triple 
interactions at z > 5, however, increases; this is due to 
the fact that the hardening timescales are now longer, as 
more massive holes create larger core radii (see eqs. 23, 
24). Overall, the number of BHs ejected in the intergalac- 
tic medium (IGM) goes up by a factor of 2. 

4.2. The mBH — er c relation 

In Figure 10 the local raBH — &c relation predicted by our 
scheme for the assembly history of SMBHs is compared to 
the fit given by Ferrarese (2002). 

The three curves in the lower corner show the raBH — & c 
relation we would obtain assuming that when halos merge 
their BHs coalesce immediately, their masses sum directly, 
and there is no gas accretion. The slope is flatter and 
the normalization lower than observed. The steepening 
of the relation for initial seed BHs in higher a peaks is 
a natural outcome of biasing. Assuming the entropy-area 
relation for the merging BHs tends instead to flatten the 
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curves. The incapacity of mergers alone to yield the ob- 
served mBH — Cc nas been pointed out by Ciotti & van Al- 
bada (2001) using an argument based on the Fundamental 
Plane. 




ct c (km s- 1 ) C7 c (km s" 1 ) 

Fig. 10. — The mBH — c c relation at z = 0. Every circle rep- 
resents one nuclear BH in a halo of given cr c . We started at z=0 
with a discrete grid of halo masses (hence, with a discrete grid of 
a c , see section 3.2 and 3.3), and performed several simulations for 
each mass. Left panel: the circles mark the results for our fidu- 
cial model with seed BHs in 3.5-cr density peaks. The solid line 
shows Ferrarese's (2002) best fit; its extrapolation to low cr c values 
is depicted by a dotted line. The holes deviating from the rela- 
tion are hosted in galaxies that experienced their last major merger 
at z > 1.5. Since then their host halos have grown due to minor 
mergers. Right panel: same for seed BHs in 4-cr peaks. The three 
curves below (from top to bottom, initial seed BHs in 3, 3.5, and 4-cr 
peaks) show the merr — o"c relation obtained assuming that when 
halos merge their BHs coalesce immediately, and there is no gas ac- 
cretion. Note how the slope is flatter and the normalization lower 
than observed. 

The scatter in our model largely reflects the time elapsed 
since the last major merger. While the final SMBH mass 
is set by the last episode of gas accretion and BH-BH coa- 
lescence, the host halo mass keeps growing through minor 
mergers and the accreting of small dark matter clumps 
with M < M res . We find that lowering the threshold for 
major mergers has the effect of smoothing the mass assem- 
bly history of SMBHs, tightening the tobh — c c relation. 
Figure 1 1 illustrates the fact that the mass-growth history 
of SMBHs does not totally reflect that of their host ha- 
los and follows a more complex pattern, with a number of 
rapid accretion episodes. 

The assembly history of two SMBHs is shown in Fig- 
ure 12; one BH (labeled as '1') ends in a main, massive 
halo at the present epoch, and the other ('2') in a satel- 
lite at z = 2.3. The final mass of the two BHs is mostly 
due to gas accretion and does not depend on the initial 
conditions, i.e., whether the seed holes are hosted in the 
3 or 3.5-cr peaks. The number of major mergers between 
halos hosting SMBHs is larger in the 3-cr peak case; only a 
fraction ends with the coalescence of the two holes before 
a triple BH interaction takes place, as gravitational sling- 
shots are now more common. As an example, for a Milky 
Way-sized halo, the fraction of major mergers that ends 
up with the coalescence of the binary without a slingshot 



drops from 75% in the 3.5-cr fiducial case to 30% for a = 3. 




redshift 



Fig. 11. — Two different realizations for the mass-assembly history 
of a galaxy halo with today velocity dispersion ctdm = 160 kms -1 
and its central SMBH. Dashed line: the halo experiences its last ma- 
jor merger at low redshift, and its hole follows the observed rrtBH — °~c 
relation. Solid line: the halo has its final major merger at z = 2.6, 
when only a small fraction of its dark mass was already in place. 
The mass of its SMBH is only about 10% of that expected from the 
m BH — °~c relation. 




redshift 

Fig. 12. — Mass-growth history of two SMBHs, one ending in a 
massive halo ('1') with <tdm = 185 kms -1 at z = 0, and one in a 
satellite ('2') with udm = 80 kms -1 at z = 2.3. Upper panel: seed 
holes in 3.5-cr peaks at z = 20. The BH mass grows after every 
major merger event due to gas accretion, independently of whether 
the other merging galaxy hosts another SMBH or not. The starred 
symbols mark the redshift of the major mergers of two parent halos 
both hosting a SMBH, but not ending up with the coalescence of 
the binary. The circles mark the redshift when two SMBHs coalesce. 
Note how most of the mass of the lighter hole is gained in the most 
recent accretion episodes. Lower panel: same but for seed BHs in 
3-cr peaks. The number of mergers not ending up with the coales- 
cence of the two SMBHs (indicated by stars) is larger compared to 
the higher bias case. The black dot marks the rare event (occurring 
in the satellite galaxy) of a decrease of the central BH mass, caused 
by the ejection of the binary in a triple BH interaction. 
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4.3. Wandering BHs 

It appears inevitable that significant numbers of triple 
BH interactions will take place at early times if the forma- 
tion route for the assembly of SMBHs goes back to the very 
first generation of stars, and even more so if the seed holes 
are more numerous and populate the low-cr peaks. As dis- 
cussed in § 3, our scheme predicts, along nuclear SMBHs 
hosted in galaxy bulges, a number of wandering BHs that 
are largely the result of minor mergers rather than of low 
energy slingshots. In practice, for minor mergers, the dy- 
namical friction timescale is longer than the Hubble time, 
and at z — the BHs are still on their way to the galactic 
center. Moreover, rare high energy slingshots can eject the 
lighter holes out of the halo, giving origin to a population 
of free-floating "intergalactic BHs" . Even rarer events oc- 
cur when the lighter BH and the binary are both ejected 
into the intergalactic medium. 

The total mass in wandering BHs ranges from 1% to 
10% of the mass of the central SMBH, for halos with 
(Tdm = 50 kms -1 and ctdm = 300 kins" 1 , respectively. 
Figure 13 shows the mass function of wandering BHs 
at z — for two typical halos, ctdm = 100 km s -1 and 
ctdm = 200 km s" 1 . 
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Fig. 14. — Contribution of nuclear, wandering and intergalactic 
BHs to the mass density in BHs as a function of redshift. Solid line: 
total mass density in all species. Dot-dashed line: nuclear BHs. 
Long-dashed line: wandering BHs retained in galaxy halos, most of 
them due to minor mergers. Short-dashed line: single and binary 
BHs ejected in the IGM after a high energy slingshot event. 
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Fig. 13. — Mass function of wandering BHs at z = 0, averaged 
over 20 Monte Carlo realizations of a galaxy-sized halo (ctdm = 
100 km s -1 , solid lines), and 20 realizations of a more massive, 
(Tdm = 200 kms -1 halo (dashed lines). The lower-left histograms 
give the contribution of slingshots to the mass function in the two 
cases considered. Errorbars are 1-cr Poissonian noise. The nuclear 
SMBHs in these halos have masses mBH = 3.8 ± 0.1 X 10 6 M Q and 
to B h = 1-0 ± 0.1 X 10 8 M , respectively. 

In Figure 14 we show the total mass in BHs predicted by 
our scheme, and the relative contribution of nuclear, wan- 
dering and intergalactic BHs, as a function of redshift. 
At z — the mass in nuclear BHs is ~ 3.5 x 10 5 M 
Mpc~ 3 , within 30% to the value given by Merritt & Fer- 
rarese (2001), somewhat larger than the value estimated 
by Yu & Tremaine (2002) from SDSS data and similar to 
an early estimate by Salucci et al. (1999). The total mass 
is dominated by nuclear BHs at every epochs, though at 
low redshifts wandering BHs become increasingly more im- 
portant. 



The cumulative local mass function of BHs is displayed 
in Figure 15. The heaviest BH hosted in a given DM halo 
has experienced several accretion episodes, lies in the nu- 
cleus, and has a well defined mass according to the tight 
mBn-a c relation. Though the total mass in BHs is domi- 
nated, in every galaxy-sized halo, by the central hole, the 
same halo hosts a multiplicity of lighter wandering BHs, 
with masses ranging from 150 M Q to approximately one 
tenth of the mass of the nuclear SMBH. Hence the to- 
tal BH mass density is dominated by wandering holes for 
masses < 10 5 M Q , and by nuclear ones above. 

Note that on group and cluster scale our scheme would 
define "wandering" all the nuclear BHs hosted in satellite 
galaxies, with only the BH of the central galaxy defined 
as "nuclear" . On the other hand, on galaxy scales, the cu- 
mulative mass contribution of wandering BHs is ~ 10% of 
the total (wandering+nuclear) mass in a given halo (Figure 
13). We extrapolate this fraction to larger scales, assum- 
ing that 10/11 of the wandering BH mass in groups and 
clusters is associated with nuclear BHs in satellite galax- 
ies, and only 1/11 is due to true wandering holes. We have 
then readjusted the mass of wandering and nuclear BHs 
in groups and clusters accordingly. 

We keep track of the positions of every single wandering 
hole within the host halos as it sinks to the center due to 
dynamical friction. The positions at the present epoch are 
typically spread between 0.1 and 0.6 of the virial radius. 
The large population of wandering BHs predicted by our 
model could be associated to the off-center, ultraluminous 
X-ray sources observed in nearby galaxies (e.g., Colbert 
& Mushotzky 1999; Makishima et al. 2000; Kaaret et al. 
2001). 

4.4. Binary SMBHs and quasars 

The fraction of halos hosting nuclear BHs is shown 
in Figure 16. When computed over all branches of 
the merger trees, i.e. considering at z = only ha- 
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los with M > 10 11 M Q with smaller and smaller ha- 
los appearing at higher rcdshifts, the occupation frac- 
tion depends on the assumed value of the resolution 
mass M ros . To avoid this problem, we compute the 
occupation fraction above a given minimum halo mass. 
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Fig. 15. — Cumulative mass density in BHs at redshift z = 0. The 
total density is dominated by the nuclear BHs contribution, though 
the contribute of wandering BHs takes over for masses < 10 Mq. 
Solid line: total mass density. Dot-dashed line: nuclear BHs. Long- 
dashed line: wandering BHs. Short-dashed line: intcrgalactic BHs. 

As discussed by Menou et al. (2001), this fraction ini- 
tially decreases with cosmic time as lower mass halos lack- 
ing nuclear BHs become more massive than the assumed 
threshold by successive mergers. Eventually, the occupa- 
tion fraction starts to increase as the total number of in- 
dividual halos drops. Note that our assumed 'BH bias' 
(seed holes in 3.5ct peaks at z — 20) assures that all ha- 
los more massive that 10 11 M Q actually host a BH at all 
epochs (triple interactions actually slightly reduce the oc- 
cupation fraction somewhat below unity). In Figure 16 
we also plot the fraction of all halos containing a nuclear 
BH that actually harbors a binary system. Summed all 
over the branches of the tree, this fraction is ~ 5 — 10% 
today. Almost all massive halos at early epochs host a 
binary, and eventually the fraction of binaries decreases 
with cosmic time because of coalescence. 

A fraction s=a 60% of SMBH binaries at z = has sepa- 
ration larger than 0.1 kpc and are still 'soft', while « 10% 
is in an advanced stage of hardening (a < 10 pc). The 
fraction of binaries strongly depends upon the total coa- 
lescence timescale, which, as discussed in § 3.4, is rather 
uncertain. We have then run a set of simulations with a 
stellar density artificially decreased by a factor of 2, which 
implies a slower hardening rate by a factor of <~ 5, on aver- 
age. In this case triple interactions are more common, due 
to the increased probability of having a binary still in the 
process of shrinking when a third intruder BH comes along. 
As a consequence, at low redshift (z < 0.1), the fraction of 
binaries becomes as large as 20%. The increased number 
of binaries is formed by close systems (a < 10 pc), hosted 
in galaxies that have not experienced any recent major 
merger. 
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Fig. 16. — Fraction of halos hosting at least a nuclear SMBH 
vs redshift (upper panel) computed weighting over all branches of 
the merger trees (dotted curve) and imposing different thresholds to 
the halo mass: M > 10 11 M (solid curve), M > 10 10 M Q (long- 
dashed curve) and M > 10 9 Mq (short-dashed curve). The last two 
quantities are not plotted at z < 0.25, when the least massive halos 
in the merger trees are already above the adopted threshold. The 
mean SMBH occupation fraction is unity at the present epoch, and, 
when weighted over all branches of the merger trees, drops below 
10 at z > 2. Black dots mark the epochs when sufficiently massive 
halos appear in the merger trees (z ~ 11 for M > 10 11 Mq, z ~ 16 
for M > 10 10 M©). Approximately 5-10% of all nuclear BHs are 
in binaries (lower panel). This fraction increases with redshift in 
massive halos. 

Binary quasars are an intrinsically rare phenomenon, as 
both SMBHs must be active at the same time. Observa- 
tionally, 16 pairs are known in a sample of <~ 10 4 QSOs, 
and among these 16 the confirmed physical associations 
are less than 10 (Kochanek, Falco, & Muhoz 1999; Mort- 
lock, Webster, & Francis 1999; Junkkarincn ct al. 2001). 
In our scheme, since only the SMBH hosted in the larger 
halo accretes and radiates during a major merger, at least 
two subsequent major mergers are required to give origin 
to a binary QSO system, the satellite BH being activated 
in the first merger event. Note that the mean timescale 
over which BHs accrete and radiate is <~ 4 x 10 7 years, 
much shorter than the typical time between two major 
mergers (Somerville et al. 2001). In our simulations we 
find a fraction of binary quasars with L > 0.01L* that is 
~ 1 — 3 x 10~ 3 at z < 4. A similar fraction is found taking a 
luminosity threshold of 0.1L* (£» was calculated at every 
redshift from the polynomial evolution given by Boyle ct 
al. 2000). If instead both BHs involved in a major merger 
were assumed to accrete, the fraction of binary quasars 
would be two orders of magnitude larger, in conflict with 
the existing data. 



5. SUMMARY 

Motivated by the recent discovery of luminous quasars 
around redshift z w 6 - suggesting a very early assembly 
epoch - and by numerical simulations of the fragmenta- 
tion of primordial molecular clouds in cold dark matter 
cosmogonies, we have assessed a model for the growth of 
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SMBHs in the nuclei of luminous galaxies out of accreting 
Pop III seed holes of intermediate masses, the cndprod- 
uct of the first generation of stars in (mini)halos collaps- 
ing at z ~ 20 from high-cr density fluctuations. As these 
pregalactic BHs become incorporated through a series of 
mergers into larger and larger halos, they sink to the cen- 
ter owing to dynamical friction, accrete a fraction of the 
gas in the merger remnant to become supermassive, form 
a binary system, and eventually coalesce. We have fol- 
lowed the merger history of dark matter halos and associ- 
ated BHs through cosmological Monte Carlo realizations 
of the merger hierarchy from early times until the present 
in a ACDM cosmology. In our scheme the current mass of 
SMBHs lurking at the center of galaxy accumulates mainly 
via gas accretion, with BH-BH mergers playing only a sec- 
ondary role. 

The main results of our investigations can be summa- 
rized as follows: 

• a simple model where quasar activity is driven by 
major mergers and SMBHs accrete at the Edding- 
ton rate a mass that scales with the fifth power 
of the circular velocity of the host halo, can repro- 
duce the observed luminosity function of optically- 
selected quasars in the redshift range 1 < z < 5; 

• hardening of BH binaries takes place efficiently both 
as a result of cuspy stellar density profiles that are 
replenished after every major merger and, to some 
extent, due to triple BH interactions. 

• although our seed BHs at z = 20 are very rare (one 
in every halo collapsing from 3.5-cr density peaks), 
the nuclear SMBH occupation fraction is of order 
unity at the present epoch. It drops to less than 
10% only at z > 2. Had we placed seed BHs in the 
4-<7 density peaks instead, the occupation fraction of 
nuclear SMBH would be approximately 0.6 today; 

• the local fraction of binary SMBHs is of order 10%, 
with half of these systems having a separation larger 



than 100 pc. (Surviving binary SMBHs have mass 
ratios 0.2 ± 0.1.) This fraction increases with red- 
shift, so that almost all massive nuclear BHs at early 
epochs are in binary systems; 

• at z < 4, binary quasars represent a fraction 1 — 3 x 
10~ 3 of all AGNs more luminous than 0.1L*; 

• the long dynamical friction timescales and BH sling- 
shots create a population of BHs wandering in galaxy 
halos and free-floating in the IGM, and contributing 
< 10% to the total BH mass density today. For a 
Milky Way-sized galaxy we estimate ~ 10 wandering 
BHs with mass between 150 and 1000 M©, and <~ 1 
wandering SMBH with 10 5 < m BH < 10 6 M (cf. 
MR). For a halo with ctdm = 200kms _1 the number 
of wandering BHs is approximately 10 times larger. 

Because of the sensitivity of our calculations on a num- 
ber of uncertain parameters, these predictions should be 
regarded only as trends. Yet, we believe our results shed 
new light on models for the assembly of SMBHs that trace 
their hierarchical build-up far up in the dark halo merger 
tree. In a subsequent paper we will explore in detail the 
possibility that the damage done to stellar cusps by binary 
BHs may be cumulative, together with the detectability of 
wandering BHs in galaxy halos, free-floating holes in the 
IGM, and coalescing SMBHs by a low-frequency gravita- 
tional wave experiment such as the planned Laser Inter- 
ferometer Space Antenna (LISA). 

We have benefitted from discussions with L. Ciotti, M. 
Colpi, D. Merritt, J. Ostriker, M. Rees, and G. Taffoni. 
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thanks the Department of Astronomy and Astrophysics, 
UCSC for kind hospitality. 
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